Synchronous solutions and their stability in nonlocally coupled 
phase oscillators with propagation delays 

Gautam C Sethi J] and Abhij it Serf] 

Institute for Plasma Research, Bhat, Gandhinagar 382 428, India 

Fatihcan M. AtajC 
Max Planck Institute for Mathematics in the Sciences, 
Inselstr. 22-26, Leipzig D- 04 103, Germany 
(Dated: May 12, 2010) 

Abstract 

We study the existence and stability of synchronous solutions in a continuum field of non- 
locally coupled identical phase oscillators with distance-dependent propagation delays. We present 
a comprehensive stability diagram in the parameter space of the system. From the numerical 
results a heuristic synchronization condition is suggested, and an analytic relation for the marginal 
stability curve is obtained. We also provide an expression in the form of a scaling relation that 
closely follows the marginal stability curve over the complete range of the non-locality parameter. 
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I. INTRODUCTION 



Phase oscillators are simple models of oscillatory activity where the essential variable 
is an angular position along a limit cycle. Systems of coupled phase oscillators provide a 
convenient paradigm for studying the collective behavior of weakly interacting limit cycle 
oscillators [lH3| and are in turn useful for modeling a wide variety of cooperative phenomena 
observed in biological, physical and chemical systems 4-sj]. One of the most well known 



models of coupled phase oscillators is the Kuramoto model [2[ which considers a system 
of globally (all-to-all) coupled set of oscillators with a spread of intrinsic frequencies. The 
oscillators synchronize to a common frequency as soon as the coupling strength goes beyond 
a threshold value. For a system of identical oscillators one achieves phase synchrony where, 
starting from an initial random distribution of phases, all the oscillators acquire the same 
phase within a time that is inversely proportional to the coupling strength. When the 
members of the set of identical oscillators have a spatial identity (e.g. when the coupling is 
of the nearest-neighbour kind or varies in strength over space) the collective states can exhibit 
a broader class of phase-locking behaviour including acquiring the same phase (synchrony) 
or a fixed amount of phase difference between neighbouring oscillators (traveling waves). 
For a nonlocally coupled continuous system of this kind it is also possible to sustain a 
more exotic state whereby phase-locked and incoherent activity can simultaneously exist at 
different spatial locations [9|, giving rise to a spatio-temporal pattern that has been termed 
chimera state [10] . 

While the original Kuramoto model is based on the idealized assumption of instantaneous 
coupling between the oscillators, there has been a growing interest in extending the model 
to include finite time delay in the coup ling; and to study the existence and stability of phase- 
locked states in such a situation llMl7|. Time delayed coupling is physically relevant in 



many real life systems and accounts for finite propagation times of signals, reaction times 
of chemicals, individual neuron firing periods in neural networks etc. One of the earliest 
investigations in this direction was made by Schuster and Wagner 11| who considered the 
collective states of a pair of delay-coupled phase oscillators and observed changes in the 
onset conditions for synchrony as well as the existence of additional (higher) branches of 



synchronized frequency states. Niebur et al. 



12] considered a two-dimensional array of phase 



oscillators that interact via time delayed nearest neighbour coupling and found that time 
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delay led to a reduction in the collective frequency of the system. Their numerical studies also 
found higher frequency states which were metastable and decayed to the lowest branch in the 
presence of noise. Kim et al [kJ studied multistability of synchronized and desynchronized 
states in a time-delayed Kuramoto system in the presence of a pinning force. More recently 
Yeung and Strogatz 15[ have carried out a systematic stability study of the synchronous 
states of a delayed Kuramoto system and have obtained analytic stability curves for the 
limiting case of coupled identical oscillators. Since the original Kuramoto model adopts a 
global coupling between the oscillators, its time-delayed version considered in all the above 
mentioned studies have assumed a single fixed amount of delay in the coupling. In a more 
realistic situation where the oscillators are spatially distributed and are coupled in a nonlocal 
fashion, such an assumption is not appropriate and one needs to consider distance-dependent 
delays. 

Models with distance-dependent delays have been proposed in the past and analyzed 
under simplifying assumptions or in convenient limits. Crook et al. [l^j considered a contin- 
uum system of coupled identical oscillators with a spatially decaying interaction kernel and 
modeled the space dependent time delay contribution through an effective phase shift term 
in the interaction. They also took the size of the system to be infinite for mathematical 
convenience but thereby effectively reduced the nature of the mutual coupling to a "local" 
one (since the interaction length is always much smaller than the system length). Zanette 
161 ] studied another simplified version of the generalized system where he adopted a distance 
independent (global) coupling between the oscillators but introduced a distance dependent 
time delay in the interaction. He obtained numerical results on the stability of synchronous 
and propagating traveling waves and also some analytic results in the limit of small de- 
lay. Ko & Ermentrout 18( recently investigated the effects of distance dependent delays in 
sparsely connected oscillator systems and found that a small fraction of connections with 
time delay can destabilize the synchronous states. We note that distance- dependent delays 



have also been considered in other related contexts, e.g., in the stabi 
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ity of equilibria of con- 



20[. To the best of our 



tinuum models of neural tissue with axonal propagation delays 
knowledge, however, there has been no systematic study of the synchronization properties 
of the full nonlocal system of phase oscillators with distance-dependent delays. Our present 
work is motivated by a desire to address this problem. The presence of both non-locality 
and time delay in the coupling pose a formidable mathematical challenge for an analytic 
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solution of the problem and one has to perforce adopt a numerical approach. In the recent 
past we have studied the generalized system to obtain chimera states in the presence of 
delay 21] . In the present work we carry out a detailed study on the existence and stability 
properties of synchronous solutions of the system. Through detailed stability curves in the 
parameter space of the mean time delay (f ) and the non-locality parameter (k) (for a fixed 
intrinsic oscillator frequency u) we delineate the relative influence of f and k on the stability 
of the synchronous state whose frequency is Q. From an extensive scan over several frequen- 
cies co, we further obtain a comprehensive stability diagram in the parameter space of Qf 
and k that is independent of oj. Our numerical results further suggest that the instability 
mechanism corresponds to a saddle-node bifurcation and the most unstable perturbation 
has the lowest mode number. These results help us obtain an analytic expression for the 
synchronization condition. The synchronization condition can be closely approximated by 
an offset exponential power law relationship between Qf and k, which is valid over a wide 
range of coupling regimes and provides a useful criterion for determining the critical value 
of the delay-induced maximum phase shift below which the synchronous state is stable. 

II. MODEL SYSTEM AND ITS SYNCHRONOUS STATES 

We consider a continuum of identical phase oscillators, arranged on a circular ring C and 
labeled by x G [-L, L] with periodic boundary conditions, whose dynamics is governed by 

• L r / mm 

dz, (1) 



^(x,t)=u-K J G(|z|)sin 



[x, t) — <j) ( x — z, t — — 

V 



where 4>(x, t) e [0, 27r) is the phase of the oscillator at location x and time t, whose intrinsic 
oscillation frequency is u > 0, K is the coupling strength and G : [-L, L] — > R is an 
even function describing the coupling kernel. The quantity v denotes the signal propagation 
speed which gives rise to a time delay of \z\/v for distance \z\ from the locations x. As 
the oscillators are located on a ring with circumference 2L, the distance between any two 
oscillators is given by the shorter of the Euclidean distance between them along the ring. In 
this configuration, the maximum distance between the coupled oscillators is L and thus the 
maximum time delay would be r m = L/v. 

We choose the coupling kernel G(|,z|) to have an exponentially decaying nature and its 



4 



normalized form is taken to be 

G(l * = 2(1.0 -e-^r "' < 2 > 

where a > is the inverse of the interaction scale length and is a measure of the nonlo- 
cality of the coupling. We make time and space dimensionless in Eqs.([T]) and ([2]) by the 
transformations t — > Kt, to — > to/K, k — > crL, z — >■ z/L, r m — > Kr m and x — > x/L and 
obtain 





di' 



(x,t) — lu — J G(\z\) sin [<f)(x, t) — (j) (x — z,t — \z\r m )] dz. 



(3) 



G(\z\) = -e- K|z| (4) 

Vl W 2(1.0 -e~ K ) V ; 

We look for synchronous solutions of Eq. that have the form: 

Mt) = + <j) . (5) 

The value of 4>q can be taken to be zero by a translation. Substituting Eq. (JSJ) in Eq. ([3]) we 

get 



CO 



G{\z\) sva.{ytT m \z\) dz 



Qr m — Qr m e K cos(f2r m ) — ne K sin(f2r„ 



(6) 



k 2 + (fir m ) 2 

which is an implicit equation in Q. Being a transcendental equation its solution can in 
principle be multi- valued in Q for a given set of parameters u, r m and k, and can lead to 



higher branches of collective frequencies as pointed out by Schuster and Wagner 
system of two coupled oscillators. 

We further define a mean delay parameter by 

G(|2;|)r m |^| dz 



llj for a 



(7) 



which weights the individual delays with the corresponding connection weights. With the 
exponential connectivity given by Eq. (jlj), this translates into 

e K — k — 1 



T = T„ 



n(e K - 1) 



(8) 



This gives values for the limiting cases: r — > for local (k — > oo) and r — > r m /2 for global 
(k — > 0) coupling. 

Fig. 1 plots the numerical solutions Q of Eq. (j^D as a function of f for k = 2.0 and several 
values of u. As the curve for u = 0.8 shows, it is possible to have multiple solutions Q for 
a given value of f . The stability of these higher states will be discussed in later sections of 
the paper. One also notes that the lowest branch shows frequency suppression as a function 
of the mean time delay f. 

We can also examine two simple limits of Eq. ([6]) which help us understand the influence 
of the non-locality parameter k on the nature of the equilibrium solutions. For k — > 0, 
corresponding to the global coupling limit, Eq. §6§ reduces to, 

n = u - 1 ~ cos( ^ rm) (9) 



This was the limiting case examined by Zanette (16J, where the transcendental nature of 
Eq. was retained, thereby permitting several multiple branches of solutions to exist, but 
the coupling was taken to be global. In the limit of Qr m <C 1, (El) gives, 

n w ^— (10) 

demonstrating frequency suppression of the collective mode due to finite time delay in this 
parameter regime. The other interesting scenario is that of k ^> 1, which corresponds to the 
situation of local coupling among the oscillators. In this case Eq. ([6]) becomes 

n = u- ■ (ii) 

Eq. (IX ip is a third order polynomial equation in Q and can at most have three real solutions 
in contrast to the higher number of multiple roots of the transcendental Eq. The local 
limit can also be approached by taking L — > oo limit in the unsealed form of Eq. ()6]) in 
which T m is replaced by L/v and k by oL. The mean delay (r^) in this system equals — 
and Eq. ([6]) can be rewritten in terms of as 

1 + iVr^ 

Crook et al. [14J considered such a limit but with an additional simplifying assumption 
regarding the contribution of the delay term: Instead of explicitly treating the delay term 
in the argument of the phase function they modeled the delay contribution by a space 
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dependent phase shift (implicit delay). Thus the basic model Eq. (JT]) was simplified in their 
case to be of the form 



^{x,t)=u-K J G(\z\) sin 



2 

(x, t) — (f)(x — z, t) H 



dz (13) 

v 

with K — 1. The corresponding equation for the equilibrium solutions is given in this case 
by 

n = "~TT^ (14) 

' oo 

which shows only a single collective frequency that experiences time delay induced reduction 
for Too < 1 and a slow rise to u for f^, > 1. 

The degree to which the non-locality parameter «; can influence the frequency of the 
synchronous equilibrium solutions of the coupled oscillator system can be seen in Fig. [2J 
This figure provides a graphical comparison of this influence by displaying plots of Q vs f 
obtained by solving Eq. with different values of k. The curve labeled NL represents 
a typical case of nonlocal coupling, here calculated for k = 2.0. The curve L is for local 
coupling with k = 10.0, whereas the global coupling curve G is calculated for k = 0.01. 
The G and L curves obtained from Eq. (Q and Eq. (TTTT) respectively are not significantly 
different from what are shown in this figure. The curve plots Eq. ( |T4|) for comparison. 
We have used u = 0.8 for the plots shown in this figure. 

III. STABILITY OF THE SYNCHRONOUS SOLUTIONS 

We now examine the stability of the synchronous solutions <pn obtained in the previous 
section. 

A. Eigenvalue equation 

The linear stability of solutions <j>o(t) of Eq. ([3]) is determined by the variational equation 
d 



-u(x,t) — — / G(\z\) cos [QT m \z\][u(x,t) — u(x — z,t — \z\r m )]dz (15) 
ot J -i 

where u(x,t) = 4>(x,t) — <frn(t). With the ansatz u(x,t) ~ e *t e vnnx ^ _\ ^ C ; ri g Z, we obtain 
the eigenvalue equation 

A = - J G(\z\) cos (nr m \z\) (l - e - A l 2 l^ e - i7m2 ) dz. (16) 



Writing A = Xr + iXi and separating Eq. (1161) into its real and imaginary parts we get 

\r = - J G(\z\) cos ( fi kkm) [1 - e" Afl|2|Tm cos (Xi\z\r m + nnz)] dz, (17) 

Xj = - J G{\z\) cos (Q\z\r m )e~ XRlzlTm sin (X I \z\T m + 7rnz)dz. (18) 

Since the perturbations u(x, t) corresponding to n = again yield synchronous solutions, the 
linear stability of the synchronous state requires that all solutions of Eq. ( Tl6l) have Xr < 
for all non-zero integer values of n. The marginal stability curve in the parameter space of 
(k, T m ) is defined by Xr = and in principle can be obtained by setting A_r ; = in Eq. (fT7|) . 
solving it for A/ and substituting it in Eq. (Tl8|) . In practice it is not possible to carry out such 
a procedure analytically for the integral Eqs. (1 1711 1 8f) and one needs to adopt a numerical 
approach, which is discussed in the next section. 



B. Numerical determination of eigenvalues 

To systematically determine the eigenvalues of Eq. ( Tl6l) in a given region of the complex 
plane we use multiple methods in a complementary fashion. First, Eq. ([6]) is solved for Q for 
a given set of values of the parameters k, lo and r m . Following that, we need to determine 
the complex zeros of the function /(A) defined as 

/(A) = A + I G{\z\) cos (nr m \z\) (l - e -M'\rm e -i^ dZj 



which is equivalent to finding solutions A of ( |T6l) . To do t 



numerical technique developed by Delves and Lyness 



lis we have primarily relied on the 



22| based on the Cauchy's argument 



principle. By this principle the number of unstable roots m of /(A) is given by 

where the closed contour C encloses a domain in the right half of the complex A plane with the 
imaginary axis forming its left boundary. Once we get a finite number for m we further trace 
the location of the roots by plotting the zero value contour lines of the real and imaginary 
parts of the function /(A) in a finite region of the complex plane (Xr, A/). The intersections 
of the two sets of contours locate all the eigenvalues of Eq. (1161) in the given region of the 
complex plane. The computations are done on a fine enough grid (typically 80 x 80) to 
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get a good resolution. A systematic scan for unstable roots is made by repeating the above 
procedure for many values of the perturbation number n and by gradually extending the 
region of the complex plane. We have made extensive use of Mathematica in obtaining the 
numerical results on the stability of the synchronous states. 



C. Results 

In Fig. [1] the solid portion of the curve shows the stable synchronous states of Eq. ([3]) for 
k = 2.0 and for various values of u and f . The terminal point on a given solid curve of Q vs 
f marks the marginal stability point. The marginal stability point is seen to shift towards 
larger values of r as one moves down to curves with lower values of u. A more compact 
representation is obtained if one plots Q — u versus Qf, since in this case the solutions 
corresponding to different values of u for a given k consolidate onto a single curve, as shown 
in Fig. [3] for k = 0.05, 2.0 and 10.0 respectively. 

It is seen from both figures that the stability domains of the synchronous solutions are 
restricted to the lowest branch where the curves are decreasing. This suggests a heuristic 
necessary (but not sufficient) condition for stability of the synchronous solutions: From Fig.[T] 
we have that dVt/df < for stable synchronous solutions, and from Eq. ([6]) we calculate 

dVl —Vtc K I 



(19) 



where 



and 



dr 1 + c K r I 

j \z\G(\z\) cos(c K Qf\z\) dz (20) 
n(e K - 1) 



Since c K > and I is bounded, the denominator in Eq. (|T9|) is positive for small values of r. 
Hence, for positive Q, the requirement dQ/df < implies the condition / > 0, that is, 

J \z\G(\z\)cxx(c K Qrf\z\)dz > (21) 

An alternative approach to arrive at the necessary condition (Eq l2~Tj) would be to make 
use of the results presented in FigFJJ We recast the dispersion relation given by Eq.flB]) in 
the form: 

Q - co = H(Qf , k) (22) 
9 



where 



H(Qf, k) = - G(\z\) sin (c K Qf\z\) dz 



It is seen from Figj3] that the stability domain of the synchronous solutions is restricted to 
the lowest branch where the curves have a negative slope. This again suggests a heuristic 
necessary condition for stability of the synchronous solutions to be H' < leading to the 
necessary condition given by Eq. (l2ip . The prime indicates a derivative of H(Qt,k) w.r.t 



As we will see later, the marginal stability curve obtained from H' or / = does lie above 
the true marginal stability curve (see Fig. EJ), confirming that condition (121j) is necessary 
but not sufficient for the stability of synchronous states. 

The points where solid and dotted lines meet in the curves of Figs. [1] and |3] mark the 
marginal stability point for the respective k values. These points are obtained for a range 
of k values and are plotted in (Qf, k) space in Fig. H] by filled points. They all lie on a single 
curve, which is analytically derived below. Our numerical results further reveal that for the 
marginal stability points, the imaginary part of the eigenvalue of the mode is zero — in other 
words the mode loses stability through a saddle-node bifurcation. It can easily be checked 
by inspection that Xj = is one of the solutions of Eq. (11 8p for any value of A#; however, it is 
not evident analytically that this is the only possible solution for Xr = 0, and our numerical 
results have helped us confirm that this is indeed the case. Hence, putting Xr = Xj = in 
Eq. (|T7|) we get the following integral relation between the parameters fl, r m and k. 



Further, we have also observed that the most unstable perturbation is the one with the 
lowest mode number, namely n — 1. Therefore Eq. (1231) with n — 1 defines the marginal 
stability curve, so the condition for synchronization takes the form 



The solid line in Fig. H] is the analytical curve of marginal stability defined by setting the left 
side of Eq. (121]) to zero, and it can be seen that the numerically calculated marginal values 
(represented by points) fit this curve perfectly. The figure also shows the stability curves 
obtained for the n = 2 and n = 3 perturbations (dashed and dotted lines, respectively) 
and these are seen to lie above the n — 1 marginal stability curve. We have carried out a 



fir. 






(24) 
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numerical check for a whole range of higher n numbers and the results are consistent with 
the above findings. 

For a system with constant delay r (i.e. if r m \z\ is replaced by r in Eq. ([3])), the cos(f2r) 
term can be taken outside of the integral in (|24p and the remaining integrand is nonnegative. 
Hence, the synchronization condition in this case becomes simply 



This agrees with the results obtained previously for constant-delay systems |15l.ll7l].Thus our 
result, as given by Eq. (123]) . generalizes the condition fl25l) to systems with space-dependent 
delays, and shows a nontrivial relation between the spatial connectivity and delays for the 
latter case. 

In order to gain some intuition into the complex interaction between connectivity and 
delays, we have obtained an approximate expression for the marginal stability curve by a 
numerical fitting procedure, yielding the relation 



for the stability of synchronous oscillations. Here, the left side involves the temporal scales of 
the dynamics (namely, it is the average time delay normalized by the oscillation period of the 
synchronized solution) while the right hand side involves the spatial scales of connectivity. 
In this view, the synchronization condition is a balance between the temporal and spatial 
scales. For high connectivity (k — > 0), the system can tolerate higher average delays in 
maintaining synchrony, and the largest allowable delays decrease roughly exponentially as 
the spatial connectivity is decreased. In the same figure we have also plotted with dotted 
curve the approximate condition fl2~TT) . which is found to lie above the marginal stability 
curve in the entire range of k. The disparity between the two curves becomes particularly 
noticeable at large values of n. 

IV. CONCLUSIONS AND DISCUSSION 

We have investigated the existence and stability of the synchronous solutions of a con- 
tinuum of nonlocally coupled phase oscillators with distance-dependent time delays. Our 
model system is a generalization of the original Kuramoto model by the inclusion of nat- 
urally occurring propagation delays. The equilibrium synchronous solutions of this system 




(25) 




fit < 0.58 + 0.56e 



-0.34k 



(26) 
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are shown to differ significantly from those of similar earlier models that had introduced 
simplifications either in the coupling or in the nature of the time delay. The equilibrium 
solutions of the lowest branch are seen to exhibit frequency suppression as a function of 
the mean time delay. We have carried out a linear stability analysis of the synchronous 
solutions and obtained a comprehensive marginal stability curve in the parameter domain 
of the system. Our numerical results show that the synchronous states become unstable 
via a saddle-node bifurcation process and the most unstable perturbation corresponds to 
an n = 1 (or kink type) spatial perturbation on the ring of oscillators. These findings al- 
low us to define an analytic relation, given by Eq. (I24p . as a condition for synchronization. 
We have also obtained approximate forms for the synchronization condition that provides 
a convenient means of assessing the stability of synchronous states. Our results indicate an 
intricate relation between synchronization and connectivity in spatially extended systems 
with time delays. 
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FIG. 1: Frequency f2 of the synchronous solutions of Eq. ([3]) for k = 2.0 as a function of the mean 
time delay f . The different curves correspond to solutions of Eq. ([6]) for different values of the 
intrinsic oscillator frequency uj. The solid portions denote stable states and the dotted ones (in 
red) unstable states. 
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2 4 6 8 10 



FIG. 2: The frequency of synchronous oscillations (G, L and NL) obtained from Eq. ([6]) for a fixed 
value of uj = 0.8 and for different values of k, representing global coupling G (k = 0.01), local 
coupling L (k = 10.0), and intermediate (non-local) coupling NL (n = 2.0). The curve -Z^oo has 
been obtained from Eq. (I14jh 
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FIG. 3: Solutions of Eq. © for the synchronous oscillation frequency for several values of k, 
plotted in terms of O — u versus Of. Note that Q — uj = H(£lf, k) for the equilibrium solutions 
(see Eq. (|22p ). In this representation, the different curves of Fig.l corresponding to different values 
of uj collapse onto a single curve for a given value of k. The solid (black) portions of the curves 
correspond to stable synchronous states and the dotted ones (red) to unstable synchronous states. 
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t 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 
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2 4 6 8 10 12 14 

K 

FIG. 4: The marginal stability curve (solid curve) in the (Clf, k) space, obtained from the lowest 
branch solutions of Eq. (I23p for n = 1. The filled circles correspond to numerical results from 
eigenvalue analysis of Eq. (|16p . and show a perfect fit to the analytical result. The dashed and 
dotted curves correspond to marginal stability curves obtained for n = 2 and n = 3 perturbations, 
respectively. The symbols S and U denote stable and unstable regions in the parameter space. 
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n — ■ — ■ — ■ — i — ■ — ■ — ■ — i — ■ — ■ — ■ — i — ■ — ■ — ■ — i — ■ — ■ — ■ — i — ■ — ■ — ■ — i — ■ — ■ — ■ — i — 1 — <-. 
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FIG. 5: A numerical fit to the marginal stability curve gives an approximate scaling law in the 
form of an offset exponential relation between £lf and k. The marginal stability curve (solid, in 
black) of Fig. [5] has been replotted along with the fitted curve (dashed, in blue). The dotted curve 
(in red) is obtained from the condition H' or / = 0. 
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